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ABSTRACT 

We consider the problem of radiative transfer in stellar atmospheres where 
the index of refraction departs from unity and is a function of density and tem- 
perature. We present modified Feautrier and Lambda-iteration methods to solve 
the equation of radiative transfer with refraction in a plane parallel atmosphere. 
These methods are general and can be used in any problem with 1-D geome- 
try where the index of refraction is a monotonically varying function of vertical 
optical depth. We present an application to very cool white dwarf atmospheres 
where the index of refraction departs significantly from unity. We investigate 
how ray curvature and total internal reflection affect the limb darkening and the 
pressure-temperature structure of the atmosphere. Refraction results in a much 
weakened limb darkening effect. We find that through the constraint of radiative 
equilibrium, total internal reflection warms the white dwarf atmosphere near the 
surface (r ^ 1). This effect may have a significant impact on studies of very cool 
white dwarf stars. 

Subject headings: radiative transfer - stars: atmospheres - stars: white dwarfs 



1. Introduction 

Problems involving radiative transfer in refractive media appear to have received lim- 
ited attention (Harris 1965; Zheleznyakov 1967; Pomranin 1968). In stellar atmospheres in 
particular, refraction is always ignored since the gas is so tenuous that the index of refraction 
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Fig. 1. — Variation of the index of refraction as a function of vertical optical depth r v in the 
nominal white dwarf atmosphere model for A = 0.948 fim. 

does not depart from unity. A notable exception is found in the atmospheres of very cool 
white dwarfs, especially those rich in helium, where the gas (or rather the fluid) density can 
reach 0.1 — 2g/cm 3 . Under these conditions, the refractive index of fluid, atomic helium can 
become as large as ~ 1.3 with very large gradients (Fig. 1) and refraction can be expected 
to affect the radiation field. In this contribution, we study radiative transfer in a stellar at- 
mosphere with refraction. While we are motivated primarily by astrophysical applications, 
the method developed here is general and can be directly applied to other problems with 
plane parallel geometry with an arbitrary, monotonic variation of the index of refraction. 

We consider refraction in a stellar atmosphere that is static, plane parallel, in local 
thermal equilibrium, in hydrostatic equilibrium, with both radiative and convective energy 
transport. The total flux is constant throughout the atmosphere and scattering is assumed 
isotropic. We further assume that the refraction of ray paths follows Snell's law (i.e. geo- 
metric optics) as is appropriate in white dwarf atmospheres. For a given run of the index of 
refraction through the atmosphere, this completely defines the trajectories of the ray paths 
and allows the derivation of a simple form of the equation of radiative transfer (ERT) in the 
presence of dispersive effects. 

There is a number of interesting physical effects that we can expect in a stellar atmo- 
sphere where the index of refraction is greater than unity. The index n u is a monotonically 
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Fig. 2. — Geometry and typical ray paths for an axially symmetric radiation field in a 
refractive, planar atmosphere with a monotonically varying refractive index n u . 

increasing function of density. In a stellar atmosphere, the index decreases from the bottom 
toward the surface where n v = 1 at r v = 0. All rays are refracted away from the upward 
vertical and some rays are internally reflected back toward the interior (Fig. 2). Intuitively, 
this should lead to an increase in the temperature profile of the atmosphere to achieve flux 
conservation. Because the largest variations of the index occur in the optically thin regions 
of the atmosphere (Fig. 1) we expect refraction to be mainly a surface effect. Furthermore, 
the group velocity of light is reduced by a factor of n v so light accelerates toward the surface 
where it reaches the speed of light in a vacuum. This directly affects the flux through its 
dependence on the velocity of propagation. The angular redistribution of light due to refrac- 
tion will reduce the limb darkening effect. Finally, the index of refraction of fluid helium is 
nearly independent of frequency from the optical to the near infrared and there should not 
be any significant chromatic effects. 

We follow the treatment of the theory of radiative transfer in dispersive media presented 
by Cox & Giuli (1968). We develop two different numerical schemes to solve the ERT in 
the presence of refraction based on the Feautrier and A-iteration methods. Both methods 
are widely used in stellar atmosphere codes and the procedures we have developed allows 
for a straightfroward implementation of the effects of refraction. To illustrate the effects 
of refraction in a stellar atmosphere, we apply both methods to solve for the radiation 
field in a realistic helium-rich white dwarf atmosphere model with T c g = 4000 K, a gravity 
of log g(cm/s 2 ) = 8 and a homogeneous composition of n(He)/n(H) = 10 6 , where n(X) 
is the number density of element X. This atmospheric structure, obtained with our white 
dwarf atmosphere code (Bergeron, Saumon & Wesemael 1995), is hereafter referred to as the 
"nominal" model. 
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In the next section, we construct the ERT for a refractive, planar atmosphere, and new 
expressions for the moments of the radiation field are given in Section 3. The numerical so- 
lution with modified Feautrier and A-iteration methods is described in Section 4. Numerical 
aspects of the solution, a detailed discussion of the effects of refraction on the radiation field, 
and the effects on the structure of the atmosphere and the emergent flux are presented in 
Section 5. 



2. The equation of radiative transfer in the presence of refractive index 

gradient 

The ERT along an arbitrary ray path is (Cox & Giuli 1968) 

= Ju{r,q) -Xuir,q)I u [r,q) + ( — 3—) (1) 

ref 



pas \ pas 



where l v is the specific intensity at point r inside the atmosphere, q is the direction of the 
curved ray path, p is the mass density, ds is an element along the curved ray path, j v is the 
total mass emission coefficient, and 

Xv = k v + <7y (2) 

is the total mass absorption coefficient, the sum of true absorption (k v ) and scattering (a u ) 
processes. The mass emission and absorption coefficients are related to their non-dispersive 
values j°, xl by 3u = and Xv = X u l n v, respectively (Harris 1965). The last term on the 
r.h.s. of Eq. (1) includes the contribution to dl v j pds coming from spatial variations of the 
refractive index. In a horizontally homogeneous medium, the geometry leads to an axially 
symmetric radiation field. In the presence of refraction, the ray path is defined by Snell's 
law: 

n u sin(9) = n' u sm{6') = constant (3) 

= 0' ± 7T (4) 

where 9 and are the polar angles coordinates with respect to the vertical z axis, parallel to 
the direction of the gradient of n v . Applying the law of energy conservation to an incident 
beam of radiation with solid angle du, propagating from a medium with refractive index n u 
into a medium with refractive index n' u , and refracted into du', 



I u cos(6)duj = r v cos{6')dw' 



(5) 
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we find that I u /nl is constant along a ray path if there are no energy losses or gains due to 
emission, absorption, or interface effects (Cox & Giuli 1968). On the basis of that assump- 
tion, the last term in equation (1) is 

dl v \ _ dl^dn^ _ ^dn^ 
ds ) re f dn u ds n v ds 

and the ERT becomes 

d ( I„(f, q) \ _ j u (r, q) - Xu(r, q)I v (r, q) ^ 



pds \ nl ) n 

In plane-parallel geometry, the derivative d/ds takes the form 

d n d d6 d , oX 

The derivative over 6 appears because the ray paths are curved. For a given ray path, 
however, Eq. (3) allows us to reduce the configuration space from {z, 9} to a one- dimensional 
curve in {z} only. The ray paths are parameterized by 8b, the value of the ray path angle at 
a reference level chosen to be the bottom of the atmosphere (Fig. 2). Hereafter, all quantities 
with subscript U B" refer to this lower boundary. A ray path is described by 

/i^j = cos Or = constant (9) 

and as for each one dimensional curve dz/ds = cos 9, we can write 

l = ™ s l < 10 > 

and rewrite the ERT (7) in a simpler form 

dV v {r v ,\iB) 

or 



I' v {t v ^b) ~ S' v {t v ) (11) 



riVBMTv)) 91 '^^ = K(t„,»b) ~ S' u (r u ), (12) 

where a primed quantity /' represents f/nl, S u = j u /Xu is the source function, o v is now 
the optical depth measured along the curved ray path, t v is the vertical optical depth, and 



/. = sign(Myi-(^f) 2 (l-^)- (13) 
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The optical depths o v and r v are related to physical distances through the relations: dr u = 
—pXudz and da u = —pXuds. 

Equation (12) reduces to the well-known expression for the ERT for a non-refractive, 
plane-parallel atmosphere if we set jj, = constant and n v — \. In the refractive case, however, 
the angle \i is now a function of /ib and t v , and at any given level, only the rays that started 
at the bottom with he given by 



will not have been reflected downward. We label the set of all possible values of /ie at a 
given t v by (. 

While Eqs. (11) and (12) are mathematically equivalent forms of the ERT, their numer- 
ical solution with finite difference schemes have very different accuracies for strongly curved 
ray paths. Equation (12) is simpler to solve numerically, as it does not require the calculation 
of a u , however we find that it is essential to integrate the optical depth along the curved ray 
path to obtain a good solution and we will use the form given by Eq. (11). 
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3. 



Moments of the radiation field 



The moments of the radiation field J u (r u ), H u (t u ), and K u {r u ) are defined as 



i 




(15) 



-i 



i 




(16) 



-i 



and 



i 
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The ERT for the non-refractive case can be rewritten in terms of these moments: 




(18) 




H v (t v ), 



(19) 



V 
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and 



&K v (t v 
Or?, 



= Jy{r u ) ~ S v {t u ). 



For the more general refractive case, the moments of the radiation field become 

1 

X f X /* c) 



(20) 



(21) 



CM 



h{hb)I' v {t v , fi B )-^-dfi B , 



and 



where 



and 



K( r ») = ^ J ^ 2 K( T u,^)dfx= - J /i 2 (/iB)^(^,/i B )^-ci/i B , 

-l 

djjb 



CM 



dfi f 



l*B 



n v ,B 



(22) 



(23) 



(24) 



(25) 



In the refractive case, H' v , and K v are more convenient quantities than J V) H v and K v 
since l' v is constant along a ray path. The construction of the moments equations (18)-(20) 
is complicated by the fact that /i and the integration domain ( are now functions of r„. We 
first consider the derivative of H' v over r„. For \i and rj given by equations (13) and (25) we 
get 



-2rj— — 



and 



d/j, 



fi 2 1 dn u 



/I n v dr v 

Taking the derivative of (22) over r„ we obtain 



(26) 
(27) 



1 d 



2 Or, 



o"^— / Vh d ^B = - ^— [r]I u ) dfi 
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, (J\~^Bmin) T > 
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a B—^Bmin 



(28) 
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Since = l' v {-HBmin) and r](fi Bmin ) = -r)(-ii Bmin ), we get 

dr v * j c 

\ J J~~( 7 ^ d ^ + \ j W~\K ~ S' v ) dfi B 



and finally 



Jl-Sl-2H' v ±-%±, (29) 



<)H " {T,) = Ut v )-S v (t v ), (30) 



which is identical to equation (18). The condition for radiative equilibrium 

B roc 



dz 



poo poo 

/ H v dv= \ X u(Ju - S u ) dv = (31) 
Jo Jo 



remains unchanged. An analogous derivation for BK' y /Br v gives 

Since r\ is antisymmetric in /i B , the last term in equation (32) is and we finally get 

BK' 1 Bn 

^jl = H' v {t v )-{ZK' v -J' v )±-^. (33) 
or v n u Ot u 

or 

^L = H v (r v )-(K v -J v )-^. (34) 
or u n u ar v 

Equation (34) contains an additional term due to refraction. In cool white dwarf atmo- 
spheres, this term can dominate near the surface, as the gradient of the index of refraction 
becomes very large (Fig. 1). This equation may be used to evaluate the radiative flux H v 
when J u , K u , and n v are known. The form of Eq. (20) that includes refraction contains a 
term in dJ y /BT v and is not useful. 



4. Solution of the equation of radiative transfer in a dispersive medium 

4.1. Feautrier solution 

In the presence of refraction, it remains advantageous to define the symmetric average 
of the specific intensity 

P'(fi B ,v,a u ) = ^[I\fi B ,u,a u ) + I'(-(j, b ,v,<t v )]. (35) 
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Differentiating with respect to a u , we obtain 

dP'{nB,v,a v ) 1 



da u 2 

and 



[I'(Hb, v, Ov) ~ I'(-fJ>B, v, <r v )\ (36) 



— = P{\iB,v,o v )-S{p v ). (37) 

The source function in a dispersive medium and under the assumption of LTE and isotropic 
scattering is (Cox & Giuli 1968) 

S' u (T v )=e u B u + {l-e v )J' u , (38) 

where e„ = k u /xu is an absorption coefficient (equivalently, 1 — e„ is the scattering coef- 
ficient, or albedo). Equations (35-37) are mathematically identical to those used of the 
non- refractive case (Mihalas 1978), and their solution will require only minor changes. 

Two boundary conditions are required to solve the second-order differential equation 
(37). At the bottom of the atmosphere, where r v >> 1, the radiation field is very close to 
thermodynamical equilibrium and 

P'b = B v (39) 

where B v is the Planck function. To ease the notation, we drop the subscript v in the 
remainder of our discussion of the solution of the ERT (Eqs. 39 - 51). 

We consider that there is no incident radiation on the top of the atmosphere (r„ = 0). 
For ray paths that exit at the surface (fi B > A*smm(0)), the surface boundary condition is 
obtained from Eqs. (35) and (36): 



dP' 



da 



= K (40) 

^surface 



Because of refraction, rays with fi B < (J> Bm in(Q) are reflected downward (Fig. 2). These ray 
paths require a different boundary condition. At the reflection point, /i — 0, I f (+(J>) = I'(—fJ,), 
and we have (36) 

dP' 

7S\ = a (41) 

a ref lection point 

We solve the equations in finite difference form, where the vertical structure of the atmo- 
sphere is divided in % — 1, 2, ... , iV horizontal layers with % — 1 corresponding to the topmost 
layer. We develop higher-order expressions for the boundary conditions using Taylor expan- 
sions: 

dP' 1 

p '» = p ' + ^r Aff « + 2' 



A<T (+ i + -(f?-S;)A<4 1 (42) 
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8P' 



da 



dP' 



&top 



(P/ - Sl)(a top - a,) (43) 



where A<jj = <7j — <7j_i and a" top is the optical depth at the boundary (either the surface or 
the reflection point) measured along the ray path. For the surface boundary condition we 
get 

= P[ + \{P[ ~ + 2a top ) (44) 

and for the boundary condition for ray paths that are totally reflected between levels r R and 
Tr-i 

= \{P'r ~ S' R )(a R + - 2a top ). (45) 

Equations (37), (39), (44), and (45) are solved by modifying the Feautrier solution 
(Mihalas 1978). The derivatives of any physical quantity / with respect to a are evaluated 
with finite differences: 

df\ ft- fi-i ^ 



dcr J ,_i crj - (7j_i 

t 2 

and for the second derivative 



(i(T 2 /i |((7 j+ i - (7j_i) 



(47) 



The ERT is solved on a discrete grid of fx B points, that defines the ray paths. Each ray path 
is labeled /iB,k, where k — 1, D. The set of radiative transfer equations (37), can then be 
written as an algebraic matrix equation (Mihalas 1978): 

- AiPU + B.P', - C,P; +1 = Li, (48) 

where Aj, Cj, and Lj are Dxl vectors corresponding to layer i, with row {A,}j,{Ci}j,{Lj}j 
corresponding to ray path j, and Bj is a -D x matrix, corresponding to layer i with elements 
{f$i}j t k=i,...,D corresponding to ray path j. The surface and bottom boundary conditions are 
constructed from equations (44) and (39). 

In the presence of refraction we need to consider the treatment of totally reflected ray 
paths. According to Eq. (45), for ray path j reflected between layers R and R — 1, the 
boundary condition for reflected rays can be written as 

£{ B *fc.*{ p '*}* - {C R h{P' R+ ih = {Mr (49) 
k 

Furthermore, the number of rays in a given layer decreases toward the top of the atmosphere. 
To maintain the dimensions of Bj, the elements corresponding to reflected ray paths are 



- 11 - 



treated differently. For angles /iB,j < l^Bminiji)-, i-e. for ray paths reflected deeper than Tj, 
we set: 

0, for j ^ k 



™>* = U io/.Lk < 50 > 

and 

{A,},- = {C ; }, = {Li}, = 0. (51) 

The solution of Eq. (47) for P'(r u , fi B ) is otherwise identical to the non-refractive case 
(Mihalas 1978): 

P; = D,P' m + v 4 

D« = (Bj — AjBj-x) 1 C, 

v, = (Bj — A i Bj_i) _1 (Lj + AjVj_i) 

D x = B^Ci 

Vl = B^Li, (52) 

and the moments of the radiation field are then obtained with 

i 

J'A T v) = \ J ' K( T v, »b) dp, (53) 



i i 



HM = ~ 2 j» ^ d/x = - J » dp, (54) 



o 



and 



i 



K(t v ) = \ J fi 2 P:(r u ,fi B )dfi. (55) 



4.2. A-iteration method 

Given the formal solution of the ERT along a ray path between optical depths av,i and 

J> w>1 ) = r^e"^ + J S' v (t)e-^) dt, (56) 
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the mean intensity (21) is given by 

J'M = \ £ />,(t-,)) d/x = A[#(f)], (57) 

where A is an operator which acts on the source function S' u . The source function (38) is 

S>„) = e„£„ + (1 - e v )A[S' v (t)]. (58) 

An iterative method can be used to find the value of the source function S' v at a given 
atmosphere level 

S' l+1 = e v B v + (1 - e v )KS[. (59) 

This solution is known as the A- iteration method (Mihalas 1978). For a discrete grid of r„ 
points, A may also be represented as a matrix operator A that acts on the source function 
vector S' (Olson & Kunasz 1987): 

J' = AS'. (60) 

The elements of the matrix A depend only on the numerical scheme adopted for in- 
tegrating over <7„ and /i to evaluate the specific intensity I' u , and the mean intensity J' u , 
respectively. A good choice is to locally interpolate the source function by a quadratic ex- 
pression in a u , which allows for analytic integration in (56) and excellent numerical precision. 
In this case the construction of the A matrix is well-known and presented in details by Olson 
& Kunasz (1987). Once the source function S' u is obtained by solving (59), V v is computed 
from (56), and the moments J' u) H' u and K' v follow by direct integration. For the A- iteration 
procedure, the boundary conditions are physically the same as in the non-refractive case: 1) 
1 1 = B u at the bottom of the atmosphere and 2) there is no incident radiation at the surface 
I' u (fx < 0) = 0. There is no need for a "surface" boundary condition for the totally reflected 
ray paths. Since the A-iteration method solves directly for I' v , we can integrate along the 
full length of the reflected ray path (Eq. 56), including the downward part. 

We found that the A-iteration method works very well for cool white dwarf atmo- 
spheres. However, it is well known that the A-iteration method converges poorly for nearly 
pure scattering media (e„ << 1) (Mihalas 1978). This difficulty can be circumvented 
with the accelerated A-iteration methods (ALI) (Cannon 1973; Olson, Auer & Buchler 1986; 
Hauschildt 1992; Hauschildt, Storzer & Baron 1993). Since the introduction of refraction 
does not affect the A-iteration method proper, the solution can be implemented just as 
easily with ALI as with the A-iteration. 
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5. Application to stellar atmospheres 



5.1. 



Numerical considerations 



As is usual, the atmosphere is divided vertically into discrete layers, r v ^ (Fig. 2) and 
the ERT is solved for a finite set of ray paths fi B ,k- Since in the refractive case the rays are 
no longer straight, the optical depth along the ray path a v must be obtained by integration 
and it is related to the vertical optical depth r v by 



where To,„ is the optical depth at the bottom of the atmosphere for upward rays and tq iV = 
for downward rays that start from the surface. We found that it is essential to calculate 
a u precisely along the curved ray paths. This is especially important in the layer where a 
given ray path is internally reflected. This is where the curvature of the ray is maximal 
and the ray can travel over a large horizontal distance (in optical depth) within that layer 
(Fig. 2). We locally interpolate //(t„) between r v grid points with a quadratic polynomial 
to evaluate a v analytically. To find the vertical optical depth r u of the reflection point 
{n{r v ) = 0), we use Eq. (13) and a quadratic interpolation of n v (r u ). The integrals for 
the moments of the radiation field are performed numerically with a quadratic (3-point) 
integration scheme. The integration consists of two parts: < fi < fi c and (j, c < /i < 1, 
because P' v (t u ,[i) is discontinuous at /i = [i c {r u ) = a/1 — n 2 v {r v = 0)/nl(T u ) (see §5.2 below). 
Rays with fi < fi c (T u ) will be reflected before they reach the surface. In Eq. (56), S' u (r u ) is 
approximated by a quadratic polynomial in each layer (Olson & Kunasz 1987). The radiative 
flux from the Feautrier solution is calculated with both Eqs. (54) and (34) to check the 
numerical consistency of the solution. 

We have verified that our methods of solutions and codes are correct in two ways. To 
check the consistency of the solutions obtained with the Feautrier and A-iteration methods, 
and their sensitivity to the resolution of the r„ and \xb grids, we solved for the radiation field 
in a white dwarf atmosphere with a fixed (T, P) structure. This nominal structure corre- 
sponds to a non-refractive helium-rich model with T eff = 4000 K, log g = 8, n(He)/n(H) = 10 6 
in radiative/convective equilibrium. In this model the convection zone extends from the bot- 
tom of the atmosphere up to r R ~ 10~ 2 . The input parameters for computing the radiation 
field are the refractive index n u (r u ) (Fig. 1), the absorption parameter e u (r u ) (Fig 3.), and 
the Planck function B v {r u ) (Fig 4.), all computed from the (T,P) structure. 

For our most precise calculation, using N = 500 layers and D = 500 ray-paths, the 
relative differences between the two methods are smaller than 0.1% in both J u and H u , and 




(61) 



-14- 



1 



r 



T 



I 1 1 1 



0.8 - 



0.6 - 



0.4 - 



0.2 - 








-6 



2 



Fig. 3. — Absorption parameter e v = k u /xv in the nominal white dwarf atmosphere model 
for A = 0.948 /im. Pure scattering corresponds to e„ = 0. 

the internal precision of both methods is better than 0.01%. The difference between the 
two solutions decreases significantly as the number of layers increases. With 50 layers, the 
moments of the radiation field can differ by as much as 1%. On the other hand, the solution 
is much less sensitive to the number of ray paths. In situations where the optical thickness 
of a layer r^j+i — r v ^ ^ 1, as is common in the deeper region of atmosphere models, the 
application of a quadratic interpolation of the source function S' v in Eq. (56) is more precise 
than the solution of Eq. (37) using a finite difference scheme. The A-iteration method is 
therefore a better choice when computing a model with a relatively small number of layers 
(N ~ 50). 

We have also compared with published solutions to a similar but simpler problem, where 
radiative transfer with refraction is considered in a dielectric slab (Abdallah & Dez 2000; 
Huang, Xia & Tan 2003). Both surfaces of the slab are maintained at fixed but different 
temperatures, radiative equilibrium is imposed, the index of refraction varies linearly between 
the two surfaces, and the absorption coefficient is constant. We note that Abdallah & Dez 
(2000) neglected the effect of dispersion on the absorption coefficient, assuming that k„ = k®, 
while the correct expression is k u = K° u /n V) where k° u is mass absorption coefficient in the 
absence of dispersive effects (Harris 1965). Setting k v = k®, we reproduce Figs. (3a-3c) 
of Abdallah & Dez (2000) perfectly. Introducing the effect of dispersion on k v raises the 
temperature profile in the slab by a few degrees. 
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Fig. 4. — Mean intensity J v (solid line), source function S v (dashed line), and Planck function 
n 2 v B v (dotted line) for non-refractive (a) and refractive (6) function of vertical 

optical depth r v in the nominal white dwarf atmosphere model for A = 0.948 /xm. 

5.2. Refraction inside white dwarf atmospheres 

We now turn to a detailed discussion of the effects of refraction in our nominal white 
dwarf atmosphere model. Figure 4 shows the mean intensity, the source function and the 
Planck function for the refractive and non-refractive cases for the same nominal (T, P) struc- 
ture.. In the refractive case, J u increases by a factor of n 2 u at large optical depth because 
Ju ~ Sv ~ n v Bv (Eq. 38). For t v < 1CT 2 , n v ~ 1 and r? v B v ~ B v . However, J„ remains 
significantly larger than in the non-refractive case because of total internal reflection. This 
can be understood from the angular distribution of the radiation. Figure 6 shows the angular 
distribution of the symmetric average of the specific intensity P' u for both cases at different 
levels in the atmosphere. Deep inside the atmosphere P' u is just the Planck function B u and 
it is the same in both cases. Toward the surface, P' v splits into two regions separated by a 
discontinuity at /i = /i c (r I/ ). The additional intensity for ji < fi c arises from the contribution 
of ray paths that have been reflected at lower r v (above the level considered). Toward the 
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Fig. 5. — Symmetric average of the specific intensity P' v = P u /n 2 v for refractive (solid line) 
and non-refractive (dotted line) cases as a function of angle \x = cos 9 at various levels in 
the nominal white dwarf atmosphere model. \x = and 1 correspond to the horizontal and 
vertical directions, respectively. The wavelength is A = 0.948 /im. 

surface, P' u doubles across the discontinuity 

P'M ~ 2P>+) (62) 

because the contribution from the source function from the optically thin regions lying 
above for /j, > /j, c is negligible but the contribution from the downward rays ([/, < 0) that 
have been reflected higher up is nearly the same as for the upward rays (// > 0), due to 
negligible absorption-emission effects inside the optically thin region. This effect should 
result in warming of an atmosphere in radiative equilibrium, because an increase in the 
mean intensity J u must generally be compensated by an increase in the Planck function B v . 
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The angular distribution of the emergent radiation, known as limb darkening, is affected 
by refraction. Limb darkening is shown in the first panel of Fig. 6 since at r = 0, Pl(fJ>) = 
P v (n) = Iu(fj)- In the non-refractive case, we have I v [t v = 0,/i = 0)/I u (0, 1) = 0.423 
for A = 0.948 /im. When refraction is introduced, the limb darkening is much weaker and 
I u (0, 0)/I u (0, 1) = 0.947. In the presence of refraction, ray paths that emerge with ji ~ 
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Fig. 6. — Radiative flux for refractive (solid line) and non-refractive (dotted line) cases 
as a function of vertical optical depth in the nominal white dwarf atmosphere model for 
A = 0.948 fj,m. The atmospheric structure is the same for both calculations. 

have been strongly deflected. The vertical optical depth r u for a ray path exiting with \i = 
is related to <r u by 



For n u (T u ) shown in Fig. 1 and o v = 1, we get t u = 0.56. This means that due to refraction, 
a horizontal viewing angle allows us to see as deep inside the atmosphere as under angle fi = 
0.56 in the non-refractive case. We can indeed see on Fig. 6 that I' v (0, 0) re f = I' u (0, 0.56) non re f. 

Despite the significant change in the specific intensity shown in Fig. 5, the radiative 
fluxes are nearly identical throughout the atmosphere (Fig. 6). This happens because in the 
absence of absorption and emission, refraction is a geometric effect that does not affect the 
value of the integral (16). By symmetry, the reflected ray paths do not contribute to the 
flux (I' u (+fi) = and Snell's law implies that nl/idpL is constant. For reflected ray 

paths in an actual atmosphere, I' u (r u ,fi > 0) > I' u (r u ,fi < 0), as emission is usually larger 
deeper inside the atmosphere and absorption processes are present. This explains why the 
radiative flux H v is slightly larger in the refractive case. For r„ > 0.1, the radiative flux 
decreases rapidly due to the presence of a convective zone, a characteristic of all cool white 




(63) 
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Fig. 7. — Pressure-temperature structure (left panel) and synthetic spectrum (righ panel) 
for the nominal white dwarf atmosphere model parameters with (solid line) and without 
(dotted line) dispersive effects. Both models are computed by imposing flux conservation 
and hydrostatic equilibrium. The filled circle indicates the level where the Rosseland mean 
optical depth r R = 1. 

dwarf atmospheres models. 

So far, our discussion of the effects of refraction on the radiation field was based on a fixed 
(non-refractive) atmospheric pressure-temperature structure. A full atmosphere calculation 
consists of finding the structure that satisfy the equations of 1) radiative transfer, 2) flux 
constancy in convection zones and radiative equilibrium in radiative zones, and 3) hydrostatic 
equilibrium. 1 Figure 7 compares such a self-consistent calculation for both the non- refractive 
(nominal structure) and the refractive cases. As expected, internal reflection leads to a 
hotter structure near the surface. The spectra resulting from these self-consistent structures 
are shown in Fig. 7. The 0.4 /im" 1 feature seen in the non-refractive case is caused by 
collision-induced absorption (CIA) by the small amount of H 2 molecules present in this 
model (n(H2)/n(He) ~4x 10 -7 ). In the refractive case, the hotter structure results in the 
dissociation of H 2 and reduces the CIA absorption. The effects of refraction become much 



lr rhe atmosphere model uses a realistic helium equation of state which predicts a very low degree of 
ionization. Electron conduction in the atmosphere is therefore negligible (Bergeron et al. 1995), in contrast 
with the models of Kapranidis (1983). 
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larger in models with lower effective temperatures or higher gravities and generally decreases 
as the hydrogen abundance increases. A detailed study of refractive effects in very cool white 
dwarfs atmospheres and spectra will be the subject of a subsequent publication. 

6. Conclusions 

Having realized that in the dense atmospheres of very cool white dwarf stars - especially 
those with a composition dominated by helium - the index of refraction departs significantly 
from unity, we have solved the problem of radiative transfer in plane parallel, LTE, static 
stellar atmospheres in the presence of refraction for arbitrary vertical variations of the opac- 
ity, temperature, and index of refraction. We show how to modify the Feautrier solution and 
the A-iteration method to solve the equation of radiative transfer in the presence of refrac- 
tion. The resulting methods are general and can be applied to any one-dimensional problem. 
For cool white dwarfs, the most important effect of refraction on the propagation of radi- 
ation in a dispersive medium is total internal reflection. This increases the mean intensity 
in the optically thin region near the surface, where the dispersive effects are strongest. If 
radiative equilibrium is imposed, this results in higher surface temperatures than in the non- 
refractive case. We also find that limb darkening is dramatically reduced by refraction. This 
first solution of the radiative transfer equation with dispersive effects in stellar atmospheres 
suggest that refraction is an important effect that must be considered when modeling the 
atmospheres of very cool white dwarfs. 
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